function a4

%  uses box method to solve
%  y'' + p(x)y' + q(x)y= f(x)   for xL < x < xR

% set boundary conditions
	xL=0; yL=2;
	xR=1; yR=1;

nx=10000;
ep=0.0001;

clf
% get(gcf)
%set(gcf,'Position', [1868 1186 573 199]);
set(gcf,'Position', [796 562 546 222]);

x=linspace(0,1,nx);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
a=zeros(1,nx-2); b=zeros(1,nx-2); c=zeros(1,nx-2);
for i=1:nx-2
	xm=x(i+1)-0.5*h;
	xM=xm+h;
	a(i)= -2 - 0.5*h*(p(xM,ep)-p(xm,ep)) + 0.25*h*h*(q(xM,ep)+q(xm,ep));
	b(i)= 1 - 0.5*h*p(xm,ep) + 0.25*h*h*q(xm,ep);
	c(i)= 1 + 0.5*h*p(xM,ep) + 0.25*h*h*q(xM,ep);
	f(i)= 0.5*h*h*(rhs(xm,ep) + rhs(xM,ep));
	end;
f(1)=f(1)-yL*b(1);
f(nx-2)=f(nx-2)-yR*c(nx-2);

% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);
y=[yL, y, yR];

plot(x,y,'-','Linewidth',1)
hold on

say=['\epsilon = ',num2str(ep)];
%text(0.03,-1.2,say,'FontSize',14,'FontWeight','bold')
%text(0.82,1.2,say,'FontSize',14,'FontWeight','bold')

box on
grid on
%axis([0 1.02 -2 2])
loc='SouthEast';
%loc='South';

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

%set(gca,'xtick',[0 1 2]);
%set(gca,'ytick',[0 1 2 3]);
%set(gca,'XTickLabel',{'0';'t_M';'2t_M'})
%set(gca,'YTickLabel',{'0';' ';' ';' '})

set(gca,'FontSize',14);
%legend(' Numerical',' Composite',' Outer','Location',loc);
set(findobj(gcf,'tag','legend'),'FontSize',14); 


function g=p(x,ep)
g=4/ep;

function g=q(x,ep)
g=13/ep;

function g=rhs(x,ep)
g=0;


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end



